close all
clear all

gavgdir = '/home/predatt/anatod/ISI/MEG/GAVG';
cd(gavgdir)
load subjname
load dataset

rootdir = '/home/predatt/anatod/ISI';
datadir = '/home/predatt/anatod/ISI/raw'; 

% input which subject to process
subjects = 1;

%% order
    % remove jumps with ft_artifact_zvalue
    % remove muscle with ft_rejectvisual: look at variance, z-scores and kurtosis
    % remove blinks & heartbeat with ICA
    % remove remaining blinks with ft_rejectvisual

%% remove jump artifacts

for sb = subjects;
infofile = strcat(subjname{sb},'_datainfo.mat'); 
subdir = fullfile(rootdir,'MEG',subjname{sb}); 
cd(subdir);
load(infofile)
cd(datadir)

fprintf('Processing subject %s... \n', subjname{sb});

cfg=[]; 
cfg = datainfo.basicinfo;
cfg.trl = datainfo.trl;
cfg.artfctdef.reject = 'complete'; % delete complete trial
cfg.continuous = 'yes';
 
% channel selection, cutoff and padding
cfg.artfctdef.zvalue.channel    = 'MEG';
cfg.artfctdef.zvalue.cutoff     = 60;
cfg.artfctdef.zvalue.trlpadding = 0;
cfg.artfctdef.zvalue.artpadding = 0;
cfg.artfctdef.zvalue.fltpadding = 0;
 
% algorithmic parameters
cfg.artfctdef.zvalue.cumulative    = 'yes';
cfg.artfctdef.zvalue.medianfilter  = 'yes';
cfg.artfctdef.zvalue.medianfiltord = 9;
cfg.artfctdef.zvalue.absdiff       = 'yes';
 
% make the process interactive
cfg.artfctdef.zvalue.interactive = 'yes';

[cfg, artifact_jump] = ft_artifact_zvalue(cfg);

% remove jump artifacts
data_no_artifacts = ft_rejectartifact(cfg);
datainfo.trl = data_no_artifacts.trl;

% find which trials had artifacts and add this information to datainfo
findsamp = artifact_jump(:,1);
remove = find(ismember(datainfo.trlorig(:,1), findsamp));
datainfo.remove.jumps = 'yes';
datainfo.remove.jumptrls = remove;

    %% save, save, save!
    
 cd(subdir)
 save(infofile,'datainfo');
 close all
end

%% remove muscle artifacts

% preprocess data, filter to see muscle artifacts
% distribute jobs across cluster, otherwise tun tmp_muscle_data function over subjects
qsubcellfun(@tmp_muscle_data, num2cell(subjects), 'memreq', 8*1024^3, 'timreq', 60*60)

% remove muscle artifacts, don't remove more than 5% or so of all trials
% check variance, kurtosis, zval

for sb = subjects

subdir = fullfile(rootdir,'MEG',subjname{sb}); 
cd(subdir);
infofile = strcat(subjname{sb},'_datainfo.mat'); 
load(infofile)
datafile = strcat(subjname{sb},'_muscle_data.mat'); 
load(datafile);
cd(datadir)
fprintf('\n\nProcessing subject %s... \n\n', subjname{sb});

    % select trials visually
    cfg = [];
    muscledata = ft_rejectvisual(cfg,muscledata);

    % remove these trials from the trl structure
    findsamp = muscledata.cfg.artifact(:,1);
    remove = find(ismember(datainfo.trl(:,1), findsamp));
    datainfo.trl(remove,:) = []; 
    datainfo.remove.muscle = 'yes';
    datainfo.remove.muscletrls = remove;

    cd(subdir)
    save(infofile,'datainfo');
    close all
    clear muscledata 

end

% delete preprocessed data files 
for sb = subjects
    subdir = fullfile(rootdir,'MEG',subjname{sb}); 
    cd(subdir);
    datafile = strcat(subjname{sb},'_muscle_data.mat');
    delete(datafile)
end        


%% perform ICA and remove heartbeat and eyeblink components

% calculate ICA components
% distribute jobs, otherwise run preproc3_ica
qsubcellfun(@preproc3_ica, num2cell(subjects), 'memreq', 8*1024^3, 'timreq', 120*60)

% choose which components to remove
for sb = subjects
    
    infofile = strcat(subjname{sb},'_datainfo.mat'); 
    subdir = fullfile(rootdir,'MEG',subjname{sb}); 
    cd(subdir);
    load(infofile)
    
    icafile = strcat(subjname{sb},'_ICA.mat'); 
    load(icafile);

    % ECG 

    close all
    % look at the timelocked/averaged components and compare them with the ECG
    figure(1); clf;
    subplot(2,1,1); plot(datainfo.ecg.timelock.time, datainfo.ecg.timelock.avg(1,:))
    subplot(2,1,2); plot(datainfo.ecg.timelock.time, datainfo.ecg.timelock.avg(2:end,:))

    % look at the spatial topography of the components
    cfg = [];
    cfg.layout = 'CTF275.lay';
    cfg.label = datainfo.label;
    cfg.comment = 'no';

    figure(2); clf;
    subplot(2,3,1); cfg.component = datainfo.ecg.ampl_i(1); ft_topoplotIC(cfg,comp)
    subplot(2,3,2); cfg.component = datainfo.ecg.ampl_i(2); ft_topoplotIC(cfg,comp)
    subplot(2,3,3); cfg.component = datainfo.ecg.ampl_i(3); ft_topoplotIC(cfg,comp)
    subplot(2,3,4); cfg.component = datainfo.ecg.ampl_i(4); ft_topoplotIC(cfg,comp)
    subplot(2,3,5); cfg.component = datainfo.ecg.ampl_i(5); ft_topoplotIC(cfg,comp)
    subplot(2,3,6); cfg.component = datainfo.ecg.ampl_i(6); ft_topoplotIC(cfg,comp)

    % EOG 

    % look at the timelocked/averaged components and compare them with the EOG
    figure(3); clf;
    subplot(2,1,1); plot(datainfo.eog.timelock.time, datainfo.eog.timelock.avg(1,:))
    subplot(2,1,2); plot(datainfo.eog.timelock.time, datainfo.eog.timelock.avg(2:end,:))
    
    % look at the spatial topography of the components
    cfg = [];
    cfg.layout = 'CTF275.lay';
    cfg.label = datainfo.label;
    cfg.comment = 'no';
    
    figure(4); clf;
    subplot(2,3,1); cfg.component = datainfo.eog.ampl_i(1); ft_topoplotIC(cfg,comp)
    subplot(2,3,2); cfg.component = datainfo.eog.ampl_i(2); ft_topoplotIC(cfg,comp)
    subplot(2,3,3); cfg.component = datainfo.eog.ampl_i(3); ft_topoplotIC(cfg,comp)
    subplot(2,3,4); cfg.component = datainfo.eog.ampl_i(4); ft_topoplotIC(cfg,comp)
    subplot(2,3,5); cfg.component = datainfo.eog.ampl_i(5); ft_topoplotIC(cfg,comp)
    subplot(2,3,6); cfg.component = datainfo.eog.ampl_i(6); ft_topoplotIC(cfg,comp)

    % look at the spatial distribution of strong ICA components (to detect eventual other artifacts)
    figure(5); clf;
    for i=1:10
        subplot(2,5,i); cfg.component = i; ft_topoplotIC(cfg,comp); 
    end;
    
      figure(6); clf;
    for i=11:20
        subplot(2,5,i-10); cfg.component = i; ft_topoplotIC(cfg,comp); 
    end;
    
       figure(7); clf;
    for i=21:30
        subplot(2,5,i-20); cfg.component = i; ft_topoplotIC(cfg,comp); 
    end;    
    
        figure(8); clf;
    for i=31:40
        subplot(2,5,i-30); cfg.component = i; ft_topoplotIC(cfg,comp); 
    end;
    
      figure(9); clf;
    for i=41:50
        subplot(2,5,i-40); cfg.component = i; ft_topoplotIC(cfg,comp); 
    end;
    
    comp_artifact = input('type the array of all artifact components: ');
    
    disp(sb)
    
    datainfo.comp.artifact = comp_artifact;
    datainfo.do_detectecgeog = 'yes';

 cd(subdir)
 save(infofile,'datainfo');

end

close all

%% check for missed blinks 
% especially important for subjects where no eye blink ICA component showed up

% preprocess, filter to see blinks
qsubcellfun(@tmp_blink_data, num2cell(subjects), 'memreq', 8*1024^3, 'timreq', 60*60)

% clean blinks manually for these subjects
for sb = subjects;
    
subdir = fullfile(rootdir,'MEG',subjname{sb}); 
cd(subdir);
infofile = strcat(subjname{sb},'_datainfo.mat'); 
load(infofile)
datafile = strcat(subjname{sb},'_blink_data.mat'); 
load(datafile);
cd(datadir)
fprintf('Processing subject %s... \n', subjname{sb});
disp(sb)

    % select trials visually
    cfg = [];
    blinkdata = ft_rejectvisual(cfg,blinkdata);

    % remove these trials from the trl structure
    findsamp = blinkdata.cfg.artifact(:,1);
    remove = find(ismember(datainfo.trl(:,1), findsamp));
    datainfo.trl(remove,:) = []; 
    datainfo.remove.blink = 'yes';
    datainfo.remove.blinktrls = remove;

    cd(subdir)
    save(infofile,'datainfo');
    close all
    clear blinkdata  
    
end

% delete preprocessed data files
for sb = subjects
    subdir = fullfile(rootdir,'MEG',subjname{sb}); 
    cd(subdir);
    datafile = strcat(subjname{sb},'_blink_data.mat');
    delete(datafile)
end        



%% the end
